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Abstract 

We statistically reconstruct a three-dimensional model of a tungsten-silver composite 
from an experimental two-dimensional image. The effective Young's modulus (E) of the 
model is computed in the temperature range 25-1060°C using a finite element method. 
The results are in good agreement with experimental data. As a test case, we have 
reconstructed the microstructure and computed the moduli of the overlapping sphere 
model. The reconstructed and overlapping sphere models are examples of bi-continuous 
(non-particulate) media. The computed moduli of the models are not generally in good 
agreement with the predictions of the self-consistent method. We have also evaluated 
three-point variational bounds on the Young's moduli of the models using the results of 
Beran, Molyneux, Milton and Phan-Thien. The measured data were close to the upper 
bound if the properties of the two phases were similar (i < E 1 /E 2 < 6). 



1 Introduction 

Predicting the macroscopic properties of composite or porous materials with random mi- 
crostructures is an important problem in a range of fields |], [l9|, fil], |38| . There now exist 
large-scale computational methods for calculating the properties of composites given a dig- 
ital representation of their microstructure; eg. permeability Q ||, conductivity || ^] and 



elastic moduli [18, 29]. A critical problem is obtaining an accurate three-dimensional (3D) 



description of this microstructure 15 



For particular materials it may be possible to simulate microstructure formation from first 
principles. Generally this relies on detailed knowledge of the physics and chemistry of the 
system, with accurate modeling of each material requiring a significant amount of research. 
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Three-dimensional models have also been directly reconstructed from samples by combining 
digitized serial sections obtained by scanning electron microscopy |26|] , or using the relatively 
new technique of x-ray microtomography |l6|]. In the absence of sophisticated experimental 
facilities, or a sufficiently detailed description of the microstructure formation (for computer 
simulation), a third alternative is to employ a statistical model of the microstructure. This 
procedure has been termed "statistical reconstruction" since the statistical properties of the 
model are matched to those of a two-dimensional (2D) image |3(], ||, ||, 34]. Statistical re- 
construction is a promising method of producing 3D models, but there remain outstanding 
theoretical questions regarding its application. First, what is the most appropriate statistical 
information (in a 2D image) for reconstructing a 3D image, and second, is this information 
sufficient to produce a useful model? In this paper we address these questions, and test the 
method against experimental data. 

Modeling a composite and numerically estimating its macroscopic properties is a complex 
procedure. This could be avoided if accurate analytical structure-property relations could 
be theoretically or empirically obtained. Many studies have focussed on this problem |fl£j l. 
In general, the results are reasonable for a particular class of composites or porous media. 
The self-consistent (or effective medium) method of Hill [22] and Budiansky [11] and its 
generalization by Christensen and Lo is one of the most common for particulate media [|1]] . 
No analogous results are available for non-particulate composites. A promising alternative to 
direct property prediction has been the development of analytical rigorous bounds (reviewed 
by Willis (47|, Hashin (T|] and Torquato Q). There is a whole hierarchy of these bounds, each 
set tighter than the next, but depending on higher and higher order correlation functions of the 
microstructure. The original Hashin and Shtrikman j2(| bounds that have been widely used by 
experimentalists implicitly depend on the two-point correlation function of the microstructure, 
although the only quantities appearing in the formulas are the individual properties of each 
phase and their volume fractions. To go beyond these bounds to higher-order, more restrictive 
(i.e., narrower) bounds, it is necessary that detailed information be known about the composite 
in the form of three-point or higher statistical correlation functions || pq |, which do appear 
explicitly in the relevant formulas. Evaluation of even the three point function is a formidable 
task, so use of these bounds has in the past been restricted to composites with spherical 
inclusions. It is now possible to evaluate the bounds for non-particulate composites |57]], 
and it is interesting to compare the results with experimental and numerical data. If the 
properties of each phase are not too dissimilar the bounds are quite restrictive and can be 
used for predictive purposes pCR . Sometimes experimental properties closely follow one or the 
other of the bounds, so that the upper or lower bound often provides a reasonable prediction 
of the actual property even when the phases have very different properties [41, 35]. It is useful 
to test this observation. 

In this study we test a generalized version [34] of Quiblier's |3(J statistical reconstruction 
procedure on a well-characterized silver-tungsten composite. Computational estimates of the 
Young's moduli are compared to experimental measurements. The composite is bi-continuous 
(both phases are macroscopically connected) and therefore has a non-particulate character. 
As such the microstructure is broadly representative of that observed in open-cell foams 
(such as aerogels), polymer blends, porous rocks, and cement-based materials. By comparing 
our computations of the moduli to the results of the self-consistent method we can test its 
utility for non-particulate media. An advantage of the reconstruction procedure we use is 
that it provides the statistical correlation functions necessary for evaluating the three-point 
bounds. Comparison of the Young's modulus to the bounds therefore allows us to determine 
the bounds' range of application for predictive purposes. 
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2 Statistical models of microstructure 



The two basic models we employ to describe two-phase composite microstructure are the over- 
lapping sphere model and the level-cut Gaussian random field (GRF) model. In this section we 
review the statistical properties of these models which are useful for reconstructing compos- 
ites. The simplest, and most common, quantities used to characterize random microstructure 
are p, the volume fraction of phase 1, s v , the surface area to total volume ratio and p^ 2 \r), the 
two-point correlation function (or 7(7-) = \p^ 2 \r) —p 2 ]/\p — p 2 ] the auto-correlation function). 
p( 2 ^(r) represents the probability that two points a distance r apart lie in phase 1. Here we 
only consider isotropic materials where p^ does not depend on direction. Also note that 
p = p( 2 )(0) and s v = -Up^{Qi)/dr. 

Realizations of the overlapping sphere model [^] are generated by randomly placing 
spheres (of radii ro) into a matrix. The correlation function of the phase exterior to the 
spheres (fraction p) is p^ 2 \r) = p'"^ for r < 2ro and p^ 2 \r) = p 2 for r > 2rQ where 

°(r) = ^(-)-±(&)\ (1) 



4 \ro / 16 \ro 

and the surface to volume ratio is s v = — 3p\np/rQ. With modification it is also possible to 
incorporate poly-dispersed and/or hollow spheres. The overlapping sphere model is the most 
well characterized of a wider class called Boolean models, which have been recently reviewed 
by Stoyan et al. 

The internal interfaces of a different class of composites can be modeled by the iso-surfaces 
(or level-cuts) of a stationary correlated Gaussian random field (GRF) y(r) (so called because 
the value of the field at randomly chosen points in space is Gaussian- distributed). Moreover, 
if r is fixed, the distribution over an ensemble will also be Gaussian. Correlations in the field 
are governed by the field-field correlation function g{r) = (y(O)y(r)) which can be specified 
subject to certain constraints [|g(?")| < g(0), lim^oo g(r) — > 0]. Invariably <?(0) is taken as 
unity. A useful general form for g is p4[ 



_ e- r /t - (r c /g)e- r/rc sin27rr/d 
9[r '~ l-(rc/0 2nr/d ' {2) 

The resulting field is characterized by a correlation length £, domain scale d and a cut-off 
scale r c . The cut-off scale is necessary to ensure 1 — g(r) ~ r 2 as r — > 0; fractal iso-surfaces 
are generated if 1 — g(r) ~ r. There are many algorithmic methods of generating random 
fields. A straight forward method is to sum (~1000) sinusoids with random phase and 
wave- vector 

[2 N 

yir) = y — ^cos(^fci-r + 0i), (3) 

where fa is a uniform deviate on [0, 2tt) and fcj is uniformly distributed on a unit sphere. The 
magnitude of the wave vectors k{ are distributed on [0, 00) with a probability (spectral) density 
P(k). The density is related to g(r) by a Fourier transform [g(r)=J^° P(k) sin kr(kr)~ 1 dk]. 
Note that P(k) > specifies an additional constraint on g(r). Although this formulation of a 
GRF is intuitive, the Fast Fourier Transform method is more efficient jl], |37f . 

Following Berk Q one can define a composite with phase 1 occupying the region in space 
where a < y(r) < (3 and phase 2 occupying the remainder. The statistics of the material are 
completely determined by the specification of the level-cut parameters and the function g(r) 
(or P(k)). The volume fraction of phase 1 is 

p = h = p/3 — p a where p 7 = (2ir)~ 2 f e~ t2 ^ 2 dt , 7 = a, (3. (4) 

J — 00 
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Table 1: The volume fraction p, two-point correlation function P2(r) and surface to volume 
ratio s v of models N, I, U and I n in terms of the properties [h, h{r) and h'(0)] of Berk's two- 
level cut Gaussian random field model. The formula h p is used for calculating the level-cut 
parameters (see Table ^) . 
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Berk and Teubner [40] have shown that the two point correlation function is p^ 2 \r) = h{i 
where 



h{r) = h' + 



., i r 9 ^ dt 



2vr7 v 7 !^ 



exp - J±- (5) 




/ a 2 - 2a[3t + /3 2 \ 
-2exp^ 2(1 _, 2) j+exp 

The auxiliary variables h and h{r) are needed below. The singularity at t = 1 can be removed 
with the substitution t = sin#. The specific surface is s v = — 4/i'(0) where 



with g given by Eqn. (|2|). 

Many more models (for which p^ 2 \r) can be simply evaluated) can be formed from the 
intersection and union sets of the overlapping sphere and level-cut GRF models. Here we 
define a few representative models which have been shown to be applicable to composite and 
porous media. A normal model (N) corresponds to Berk's formulation. Models can also be 
formed from the intersection (I) and union (U) of two statistically identical level-cut GRF's. 
Another model, l n , formed from the intersection of n primary models, has also been found 
useful. The statistical properties (p, s v and p^) of each model are given in terms of the 
properties of Berk's model [Eqns. (|), © and (§)] in Table | @. 

Since the volume fraction of the models is a function of both level-cut parameters (a,j3) 
there is a continuum of choices which correspond to a given volume fraction. For example 
( a ,/3)=(-oo,0.84), (-0.84,-0.25), (-0.25,0.25), (0.25,0.84) and (0.84, oo) in Berk's model 
all correspond to p = 20%. The final two choices are statistically identical to the first two and 
therefore provide nothing new. We note that a small change in these parameters will only 
slightly alter the microstructure, so as a compromise between simplicity and generality it is 
suggested only three distinct cases be considered: (i) the common single-cut field (a = — oo); 

(ii) a symmetric two-cut field (a = —0) and (iii) an asymmetric two-cut field. A concise way 
of expressing this is to take p a = | — %h p and pp = | + (1 — where h p for models N, 
I, U and I n is given in Table |l] and c G [0, 1]. Setting c=0, 1 and ^ gives cases (i), (ii) and 

(iii) respectively. The implicit formula for finding (a, (5) from p a and pp is shown in Eqn. 

As an example the results of the calculation for nine different models (N, I and U at c=0, ^ 
and 1) at volume fraction 20% is shown in Table ||. The model N (c=0) is the single level 
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Table 2: The level-cut parameters a and (3 for different Gaussian random field models are 
calculated by solving Eqn. (||) where p a = | — %h p and p ( g = | + (l — %)h p . h p is shown in 
Table |l| for each model and c = 0, i, 1. This table shows the results of the calculation for 
volume fraction p=20%. 
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c= 1/2 
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Normal (N) 


-oo -0.84 
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-0.25 0.25 


Intersection (I) 


-oo -0.13 


-1.09 0.22 
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-oo -1.25 
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cut GRF previously used by Quiblier pH] and Teubner pQ]. Model I (c=l) has been used to 



model aerogels [33] and model N (c=l) is Berk's M two level cut model of microemulsions. 



3 Statistical reconstruction 

3.1 The Joshi-Quiblier-Adler (JQA) Approach 

The basic idea of statistical reconstruction is to generate a three dimensional (3D) model of a 
random composite using only statistical information measured from a two dimensional (2D) 
image. Since 2D cross-sections are readily obtained by common experimental methods this 
provides a very attractive method of modeling porous and composite media. Quiblier pC| ] 
developed a method capable of producing a three-dimensional model with a specified volume 
fraction and two-point correlation function. The method was first studied in two dimensions 
by Joshi [^J] and has been extended by Adler Thus it has been called the JQA model. We 
first give a summary of the procedure to demonstrate its equivalence to the single-cut GRF 

(2) 

model discussed in the last section. First p and p e ^pt are measured from an experimental 
image. The level-cut parameter (3 is fixed by the volume fraction [Eqn. with a = — oo]. 

(2) 

Second (fexpt^i) is obtained at the set of discrete points where Pgxpt is measured by inverting 

Eqn. (||) with the left hand side set to Pexpt( r i) an d ot = — oo. To carry out the inversion the 
right hand side of Eqn. (||) is expanded as a series in powers of g; 

oo 

p& in)=P 2 +p{i-p)Y. c l kfcxpt {n)] m (7) 

m=0 



where the coefficients C m depend on certain integrals of Hermite polynomials pOJ . The series 
converges very slowly for g near one, and Alder [|I| has discussed the inversion procedure 
in this case, and possible situations where the equation has no solution. We note that the 
inversion can also be simply carried out by numerical integration of Eqn. (|B|) and a standard 
non-linear equation solver. 

Next it is necessary to generate a GRF with field-field correlation function g e xpt(?")- 
Quibler's original formulation involved the solution of a very large system of non-linear equa- 
tions for the terms of a convolution operator used in his definition of a GRF. Adler 
simplified the procedure by reformulating the problem in terms of Fourier transforms, which 
is equivalent to a three-dimensional version |37]] of Rice's |3^] method for Gaussian processes 
in one dimension. In terms of the definition given in Eqn. (Q) the inversion is equivalent to a 
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numerical integration of P(k) = ~ g(r)r k sin krdr, where g is known only at a discrete set 
of points. The inversion methods described above for g(r) do not guarantee that P(k) (or its 
equivalent in other formulations) is greater than zero. However Adler Q found that if P(k) 
is negative at some points it is also small, and can be replaced with zero. Finally the GRF 
is thresholded in the usual way to obtain the reconstructed microstructure. Thus the JQA 
method produces a single-cut Gaussian random field, which we have termed model N(c=0). 

In this paper we employ a different implementation of the JQA method [p4fj . At this stage 
we restrict attention to model N(c=0). First, the volume fraction of the model is set to that 
of an image: p m od=Pexpt- Second, the experimental two point correlation function is fitted by 
varying the morphological parameters of a given g(r) [Eqn. (0)] to minimize the non-linear 
least squares error 

V f \rP (r ) - r? I 2 
2^i=lLPexpU'«7 PexptJ 

where Nj is the number of experimental points to be fitted. Numerical integration is used 

to find Pmod( r «) [Eqn. @]. The minimization is very fast, but several starting points should 
be used as a check against local minima. Once the parameters of g(r) are known an analytic 
form of P(k) is used to generate the coefficients of the GRF (||). The reconstructed model 

(2) 

is obtained by thresholding the random field as described above. Note that p m od( r i) win n °t 

(2) 

match Pex P t( r «) exactly at each point as would be the case with Quiblier's procedure. However 
with this choice of g(r), P{k) is guaranteed to be positive. More general functional forms of 
g(r) can also be employed. 

For isotropic materials p and p^ (r) can be exactly measured from a two (or one) dimen- 
sional image. Hence the application of the JQA method results in a model which shares p, s v 
and p^ 2 ) (r) with a real composite. The question is whether or not the model provides an accu- 
rate and useful representation of the original microstructure. In certain cases it appears to, in 
other it does not. First, predictions of transport properties (conductivity and permeability) 
obtained from reconstructed porous models under-estimate experimental and numerical data. 
Second, the percolation threshold (the volume fraction at which the pore space or inclusion 
phase is no longer macroscopically connected) of model N(c=0) is around 10% [[37]] for the 
model, but many materials exhibit lower thresholds fl36f| . Both points indicate that the pores 
(or inclusions) of the model are not sufficiently well connected to mimic many physical ma- 
terials. Third, we can test the model by trying to reconstruct several of the distinct models 
defined in the previous section. The results are shown in Fig |l[ Even though model N(c=0) 
is able to reproduce the correlation functions reasonably well, the reconstructions (i) do not 
in all cases appear to reproduce the original microstructure and (ii) look quite similar to one 
another. This indicates that irrespective of g(r), and the original image, model N(c=0) can 
only generate microstructures that are similar to those shown in the final row of Fig. ffl. 

Quiblier [ 30 1 has suggested that the ability of the method to generate a reasonable model 
depends on the validity of the hypothesis that: all the necessary information about the mor- 
phology is contained in the auto- correlation function. Suppose this were true. Then since the 
JQA method is sufficiently general to re-produce all reasonable two-point correlation functions 
(see Fig. Ill), it must also be able to generate all types of morphology. The discussion above 
(and Fig. |l|) indicates that model N(c=0) can only produce a limited class of microstructure 
and therefore that the hypothesis is false. This does not mean that the the method cannot 
produce useful models, but we argue it will do so only when the original material is approx- 
imately contained in the same limited class. If it is not, then a model from a different class 
needs to be considered. We discuss this issue in the following section. 
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Figure 1: The statistical reconstructions (bottom row) of four different two-dimensional im- 
ages by adjusting the parameters (r c , £, ci and /3) of model N(c=0) to fit the auto-correlation 
functions (middle row) of the original models (shown in the top row). The procedure used 
is very similar to that of Quiblier. The results suggest that model N(c=0) cannot mimic all 
types of microstructure, even though it can reproduce the two-point correlation functions of 
all the cases shown. The models in the top row are (from left to right), overlapping spheres, 
models N(c=l), I(c=l) and U(c=l) (see Sec. |2|). The images are 128x128 pixels, and the 
length scale of the correlation functions shown in the 2nd row extends to 32 pixels. 



3.2 Using several models 

From the foregoing discussion it is clear that a model more general than a single level-cut 
random field is needed to reproduce the random isotropic microstructures seen in many com- 
posite materials. At present there is no one model that can achieve this. Instead it has been 
proposed that a number of morphologically distinct models (N, I and U of Sec. ||) be incor- 
porated |3"3|| . This is very simple to do by using the relevant formula for p^ od {r) (Table ||) 
in Eqn. (^|). It was found that many of these models were able to match any given Pexpt( r ) 
(providing further evidence that the two-point correlation function does not provide sufficient 
information for a useful reconstruction). The problem then becomes how to choose the best 
model. Clearly higher-order statistical properties need to be taken into account. 

The quantities p and p^ (r) are the first and second of an infinite hierarchy of correlation 
functions, the three-point function p( 3 \r, s,t) representing the probability that three points, 
distances r, s and t apart fall in phase 1, and so forth. Two random composites can only be 
said to be statistically identical if their iV-th order correlation functions (N = 1, 2, 3, 4 . . .) are 
identical |48|]. Therefore the exact statistical reconstruction of a composite requires matching 
all of the correlation functions^} An obvious method of choosing the best of several models 
is then to compare of each model to experimental data. As well as being memory and 
time intensive || it has been shown [34] that p^ 3 \ like p^ 2 \ may not contain the relevant 



1 Technically, the model and composite may still differ by a point-process [^9[, but this is unlikely to effect 
the macroscopic properties. 
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morphological information (ie. it does not provide a strong signature of microstructure). A 
comparative study of several high-order statistical quantities found that the simplest and 
most discriminating signature of microstructure was the chord-distribution functions for each 
phase pW(r) (j = 1,2). pW\r) is the probability that a randomly chosen chord in phase j 
has length r. A chord is defined as any line-segment which lies entirely in phase j with end 
points at the phase interface. Like p^ the chord functions are the same whether measured 
from a two or three dimensional element of the microstructure. 

The chord-functions can be employed in a reconstruction algorithm as follows. First the 
morphological parameters [the length scales in Eqns. (|l|) and (pi)] of each model (overlapping 

(2) 

spheres, or the the various level-cut GRF's) are chosen to fit p e J pt . We have found that most 

models are able to provide a reasonable fit of Pcxpt( r ) ( e -§- Ep^ < 0.1). If this is not the case 

the model is unlikely to provide a useful reconstruction and may be rejected. Second, of the 

(i) 

remaining candidates, the model that best reproduces the experimental chord functions Pexpt 
is selected as the best reconstruction. We quantify the error by a normalised least square 
sum; 

where /?rec are the measured chord-distributions of the reconstruction for phase j = 1,2 (at 
M points). The final reconstruction thus has approximately the same chord functions as the 
experimental image as well as sharing the low-order quantities p and p( 2 \r). 

A limitation of this 'model-based' technique is that one of the of models tested must, 
for some choice of its morphological parameters, be able to approximately reproduce the 
experimental microstructure. For example, it would be unlikely that a model derived from 
the iso-surfaces of a random field would be able to mimic the highly structured morphology 
of randomly packed hard spheres. In such a case we would expect that none of the model 
chord functions would reproduce the experimental data. At present there have not been a 
sufficient number of studies to provide numerical criteria on EpV' for acceptability. The 
general approach we have outlined is not restricted to the models given in Sec. ||. Ultimately 
it would be useful to incorporate poly-disperse overlapping spheres and other Boolean models 
such as those based on Poisson and Voronoi polyhedra. The latter models have proved useful 



in the analysis of miner alogical materials [25] and flow in porous filters [10]. There may also 
prove to be more useful discriminants of microstructure than the chord-functions. 

To conclude this section we note that a recent study [Q has considered a 'model inde- 
pendent' scheme based on sequentially moving filled pixels (representing the target phase) 
on a grid so that the reconstruction reproduces statistical properties of the original im- 

(2) (2) 

age. The method differs from ours in that numerical estimates of pr ec replace p^od m 
Eqn. (||), and Eqns. @ and (||) are coupled. The authors also employ the lineal path 
distribution function L^\r) in Eqn. (g) which is related to the pore-chord functions by 
p U)(r) = \s v d 2 L^{r)/dr 2 g|. 

4 Elastic properties 
4.1 Theory 

The basic information required for the evaluation of the effective moduli are the volume 
fractions, and elastic moduli, of each phase; p (fraction of phase 1), q = 1 — p, and pi 
(i = 1,2). A common approximation for the effective moduli is the self-consistent method 



(SCM) of Hill [^] and Budiansky [11] which involves solving the equations of elasticity for 



S 



a spherical particle of phase 1 surrounded by a medium of unknown effective moduli K e fi e . 
The results are obtained by solving 



-?— + —*— = (10) 

K e — K 2 H e — oK e + 4/i e 

P | Q = 6(K e + 2/i e ) 



He - [12 He~ Hi 5yLt e (3K e + 4^ e 

for K e , He- in the case where one of the phases is perfectly soft or rigid the results ex- 
hibit a percolation threshold of p = i. The formula is also symmetric to phase interchange 
[n e (ki,gi, &2,<72,p) = K e (k2, 52) k\, gi, 1 — p) etc.]. These facts limit the applicability of the 
SCM, since most composites have lower percolation thresholds and many are not symmet- 
ric Jl9|. A more realistic formula is obtained using a generalized SCM (GSCM) g| for 



the case of a particle of phase 1 surrounded by a spherical shell of phase (embedded in a 



medium of the effective moduli). The result is complicated [12] and not reproduced here. 
The GSCM has zero percolation threshold, and is not symmetric under phase interchange. 
For non-particulate media it is not clear which phase should be associated with the inclusions 



and which with the matrix. Below we consider both cases. Christensen [12] found that the 
GSCM provided a better prediction of composite properties than other common methods, so 
we shall not consider these here. It should also be noted that the volume fraction is the only 
microstructural information included in the SCM and GSCM results. This means that these 
formulae are insensitive to the distribution of each phase, or rather that each formula has a 
"built-in" microstructure, which may or may not match the experimental one. 

The difficulty in deriving general theoretical results for predicting the elastic properties 
of random composites has provided the impetus for the development of rigorous bounds 0]. 
For orientationly isotropic materials the bounds take the general form [|31] J, 

<o - lJ ^S^X \ " l < «e < w - 3p ^;:f <i2) 




(M-')- ™,'^' ) <*..<W- 6W 6 ( ( ^g ) (13) 

where for a variable b, {b} = pb\ + qb 2 and (b) = qb\ +pb2- The additional parameters T, A, S 
and depend on the level of microstructural information available. If any of the moduli (nt 
or Hi) are zero then the lower bounds vanishes. Similarly if any of the moduli are infinite the 
upper bound diverges. 

If only the volume fractions of the composite are known 

T = ^\ A = M2 , 5= ; i + 2/Xl , & e = M2(9«2 + 8^ 2 ) (14) 

/Ui(9ki + 8/ii) k 2 + 2/i 2 



and Eqns. ( |12|) & ( p~3|) are the bounds of Hashin and Shtrikman [£0| for the case Hi — Hi 
and K2 > hi- The bounds only apply to well-ordered materials [(k,2 — ki)(h2 ~ Hi) > 0] 
and the inequality signs in the bounds must be reversed if H2 < Hi an d K 2 < «i- If further 
information is available in the form of three-point statistical correlations it is possible to derive 
more restrictive bounds in terms of the microstructure parameters 



9 f°°dr f°°ds f 1 An D u s ( (3) n ^ p (2 \r)p^ 2 \s) 



Ci = — - - / duP 2 {u) pW(r, a ,t)- £ v ^ w (15) 
2pq Jo r Jo s 7-1 \ P J 

5 , 150 f°°dr f°°ds f 1 , „ . , / m . . p^(r)p^(s)\ , . 

I I I ' a (tt) p (3) (r,M) - W W ■ (16) 



5 , 150 (°°dr f°°ds f 1 , „ , 

m = tttCi + / — / — / dnP 4 ( 
21 7pq Jo r Jq s V-i 
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where t 2 = r 2 + s 2 — 2rsu and P%{u) = ^(3u 2 — 1) andP^u) = |(35w 4 — 30u 2 + 3) are Legendre 
polynomials. In this case we have 

r = (M ) CJ a = (m)c5 e= {2k ., )( + 5{ , )v 

(128k- 1 +99/x- 1 ) c + 45(//- 1 )^ ' 1 j 

Here we have used the standard notation (b)^ = Ci&i + C2&2 and (6)^ = 77161 + 77262 where 
6j is any function of /ij and Kj, C2 = 1 ~~ Ci an d % = 1 — Vi- With the parameters given in 
Eqn. (|l7|), the bounds on k are due to Beran and Molyneux @ while the bounds on /i are 
those of Milton and Phan-Thien [p9]. The development of these bounds have been recently 
reviewed |4l], |3l]]. Below we consider the Young's modulus [E = 9k/// (3k + /j)] and Poisson's 
ratio [v = (3k — 2/i)/(6« + 2/i)] of a composite. The bounds on E |l^] and 1/ [|5(| are 

9ki(j,i 9K u fi u 3ki - 2/i n 3k m - 2/j/ 

< E e < - & — < v e < — , (18) 



3k ; + Hi 3k u + /i u 6ki + 2/i n 6k„ + 2/j z 



where the subscripts refer to the upper and lower bound of Eqns. ( |T2"D &: (|T^). Depending on 
the level of microstructural information employed to find the bounds on k and /j we refer to 
Eqn. ( |l8| ) as the Hashin and Shtrikman (HS) or Beran, Molyneux, Milton and Phan-Thien 
(BMMP) bounds respectively. The microstructure parameters £ and 77 have been evaluated 
for hard and overlapping spheres |^] , level cut Gaussian random-field models [ 37 1 and can be 
evaluated for other Boolean models |j23|| . 

4.2 Computation 

A microstructure made up of a digital image is already naturally discretized and so lends 
itself to numerical computation of many quantities. For computing elastic moduli, there are 



two methods available: a finite element method [18], and a finite difference method [29]. The 
finite element method uses a variational formulation of the linear elastic equations, and finds 
the solution by minimizing the elastic energy via a fast conjugate gradient method. The 
finite difference method formulates the linear elastic equations directly in a finite difference 
approach, and solves the resulting set of linear equations with a similar conjugate gradient 
method. 

For a porous material, with one solid phase and one pore phase, either method can be 
used, as the zero normal force boundary condition at a solid-pore boundary is easy to handle 
in either method. When there are solid-solid boundaries between two different phases, the 
boundary conditions become continuity of displacement and continuity of normal force. This 
is harder to implement in the finite difference method, while it is just as easy to do as in 
the solid-pore case for the finite element method. We have used the finite element method 
exclusively in this paper. 

The finite element method is one that has been especially adapted for digital images. It 
is for linear elasticity only. Each pixel, in 3-D, is taken to be a tri-linear finite element []14|]. 
There can be any number of phases, whether isotropic or anisotropic, as long as each phase can 
be adequately depicted within the resolution of the digital image used, and can be described 
with a single elastic moduli tensor. Thermal strains can also be easily handled [fT^I . The 
digital image is assumed to have periodic boundary conditions. A strain is applied, with 



the average stress or the average elastic energy giving the effective elastic moduli [41, Q 



Details of the theory and programs used are reported in t he papers of Garboczi fc Day |jj 



and Garboczi [17]]. The actual programs are available at http://ciks.cbt.nist.gov/garboczi/ 
Chapter 2 
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Figure 2: Cross-sections of (a) the overlapping sphere model, and (b) the best reconstruction 
(model Iio c=0). The volume fraction is p=20% and the images are 96x96 pixels. 



4.3 Test case: Overlapping sphere model 

Before going ahead and reconstructing the W-Ag composite we first tested the procedure for 
the over-lapping sphere model for which the statistical properties can be analytically evalu- 
ated [44|. A 3-D realization of this model is also easy to generate, so that the reconstruction 
can be carefully tested. The overlapping sphere model [Fig §(a)] has previously been recon- 
structed using the models derived from the level cut GRF's |54| . Model Iio (c=0) was found 
to be the best reconstruction [Fig §(b)]. To gauge the accuracy of the procedure for the elas- 
tic properties of the tungsten-silver composite (see the following section) we have computed 
the elastic moduli of the overlapping sphere model and its reconstruction at volume fraction 
p = 20% (of the phase outside the spheres). The moduli of each phase at each temperature 
is set to the corresponding value for silver and tungsten. The results, shown in Fig. ||, in- 
dicate that the procedure performs very well. When the silver has non-zero elastic moduli, 
for temperatures below the melting point of silver, the reconstruction provides an extremely 
good prediction of E (error <1%). At temperatures above the melting point of silver, the 
silver is taken to have zero shear and bulk moduli, and the reconstructed model is 9% stiffer 
than overlapping spheres. Since the moduli depend most strongly on microstructure at high 
contrast (contrast = the ratio of the Young's moduli between the two phases) the latter error 
is likely to be more indicative of the ability of the model to reproduce the microstructure of 
overlapping spheres. Nevertheless the reconstruction provides a reasonable model. Similar 
agreement was seen for the Poisson ratios (not shown). 

An advantage of the method is that the BMMP bounds [Eqn. (El|)] can be evaluated 

(3) 

(since p^ec can be computed). The microstructure parameters for the reconstruction [model 
Iio ( c= 0)] were estimated as £i = 0.43, rji = 0.35 Q]. The bounds (see Fig. ||) are significantly 
more restrictive than those of Hashin and Shtrikman, and are seen to bound the finite element 
moduli of both the overlapping sphere and reconstructed models. Above the melting point 
the lower bounds are identically zero, as we have taken the elastic properties of silver to be 
zero at this point. Actually, above the melting point of silver, it would be more reasonable to 
suppose that the silver has a non-zero bulk modulus, with a zero shear modulus. This is only 
important when comparing with experimental results, however, and not in this model-model 
comparison. 



5 Application to a tungsten-silver composite 

To quantitatively test the reconstruction method, experimental data need to be available 
giving a picture of the material, properties for each individual phase, and overall composite 
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Figure 3: Young's modulus of the overlapping sphere model (OS) compared with data obtained 
from the best reconstruction (Recon.) [model Iio (c=0)] (Finite element method). The Beran, 
Molyneux, Milton and Phan-Thien (BMMP) bounds are seen to be more restrictive than the 
Hashin and Shtrikman (HS) bounds for both models. 



properties. A well characterized system, suitable to test the reconstruction procedures, is 
provided by the tungsten-silver (W-Ag) composite of Umekawa et al. [fjq| . This composite was 
produced by infiltrating a porous tungsten solid with molten silver (volume fraction of silver 
p=20%). The Young's modulus of the composite was measured at a range of temperatures 
above and below the melting point of silver (960°C). The elastic moduli of each phase were 
obtained by measuring the moduli of pure samples of tungsten and silver at each temperature. 
This data cannot be used directly because both phases of the composite actually contained tiny 
spherical pores. These will reduce the Young's moduli and Poisson's ratio of each phase. This 
effect can be accounted for by applying well known results for dilute spherical inclusions [jE| . 
For porous materials (porosity 1) the formulae can be rewritten as 

^ = E m -<f>E m ( 9 -* Um - 5u2 A (19) 
y 7 - 5i/ m J 

U * ~ Um -2^{ 7^ ) (20) 

where E m , u m denote matrix properties and Eq, are the porosity modified values. The 
tungsten matrix had an internal porosity of 1% while the silver phase had a porosity of 10% 
at room temperature, decreasing linearly to 5% at the melting point. Table [|| shows the phase 
moduli used at different temperatures. 

To reconstruct the W-Ag composite we digitize a photograph of the sample (Fig. ||). 
All points below a selected threshold grey-level are set to black (the silver phase) while the 
remainder is set to white. The image was blurred and re-thresholded to remove the pores of 

(2) 

the tungsten matrix (which appeared as silver) . This had little effect on p e xpt and Pexpt > but 
a significant effect on the measured silver chord distribution. The resulting image is shown 
in Fig. ||. The image actually has a silver content of only 13.5%, significantly lower than the 
nominal value of 20%. The statistical properties of the image are compared with those of 
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Table 3: The moduli of the silver and tungsten phases, as a function of temperature, after 
being corrected for the internal porosity of each phase. 





Silver 


Tung 


sten 


Temp(°C) 


E(GPa) 


V 


E(GPa) 


V 


25 


71 


0.36 


400 


0.28 


200 


69 


0.36 


392 


0.28 


400 


63 


0.36 


383 


0.28 


600 


54 


0.36 


373 


0.28 


800 


45 


0.37 


363 


0.28 


860 


42 


0.37 


361 


0.28 


910 


39 


0.37 


359 


0.28 


950 


37 


0.37 


357 


0.28 


960 


37 


0.37 


356 


0.28 


960 





0.50 


356 


0.28 


1020 





0.50 


354 


0.28 




11 trial reconstructions in Table while 2-D slices of the models themselves are shown, for 
purposes of illustration, in Fig. |6[ Several of the models were unable to reproduce Pexpt anci 
were considered no further. A comparison of the chord distributions indicated that model 
N (c=0) provided the best reconstruction. The auto-correlation function of this model is 
compared with experimental data in Fig. [?]. The experimental and model chord-distributions 
are shown in Fig. ||. The silver distribution Pexpt i s wen reproduced by the model at all lengths 

(2) 

shown, while p^od performs less well at small chord lengths. Two and three dimensional 



images of the model are shown in Fig. M (shown at the same scale as Fig. and in Fig. 10 



For the purposes of computing the elastic properties of the model we maintain the length 
scale parameters (£, r c and d) and alter the level cut parameters (a and (3) of the model such 
that p rec =20% (in accord with the experimental composite). The Young's modulus, computed 
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Figure 5: The cropped (640x640 pixels) binary image obtained from the scanned image. This 
is the sample used to calculate the statistics of the composite (side length 198. 7^m). The 
black phase corresponds to silver. 



Table 4: A comparison of the statistical properties of 11 reconstructions with those of the 
experimental composite; p=13.5% and s v =0.17fim~ 1 (obtained from Fig. ||). Many of the 
models are able to reproduce the low order statistical properties of the composite (Ep^ < 
0.1). This shows that p( 2 \r) does not uniquely specify composite microstructure. 



Mod. 


c 


r c 


e 


d 


■S v 








N 





2.16 


2.15 


13.0 


0.19 


0.05 


0.03 


0.28 


N 


1 

2 


28.1 


28.2 


22.0 


0.18 


0.11 






N 


1 


CO 


CO 


25.6 


0.18 


0.10 






I 





2.88 


2.89 


12.5 


0.20 


0.05 


0.28 


0.88 


I 


1 

2 


13.4 


13.5 


15.9 


0.18 


0.06 


0.32 


0.53 


I 


1 


32.7 


32.1 


17.4 


0.17 


0.05 


0.11 


0.38 


U 





2.69 


2.68 


13.1 


0.18 


0.05 


0.08 


0.30 


u 


1 

2 


60.2 


144 


35.5 


0.20 


0.15 






u 


1 


CO 


CO 


43.0 


0.20 


0.15 






Iio 





4.75 


4.76 


12.6 


0.21 


0.06 


0.33 


0.58 


OS 






r =3.75 




0.21 


0.25 


0.49 


0.45 



using the finite element method, is compared with the experimental data of Umekawa et al. 



in Fig. 11. For the temperature region below the melting point of silver the maximum error 
is 4%, a very good result. Above the melting point of silver, when the silver phase is taken 
to have a zero bulk and shear modulus, the error is only 3%. The agreement may actually 
be better than that, however. Since the elastic measurements were dynamic measurements, 
the liquid silver can be considered as being trapped on the time scale of the experimental 
measurement, before any significant flow could take place, and so could still contribute to the 
effective moduli via its non-zero liquid bulk modulus. Just before melting, the silver had a 
bulk modulus of about 35GPa. If we take the bulk modulus to be somewhat lower, in analogy 
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Figure 6: Cross-sections of a portion of the original image (a) and the eleven trial reconstruc- 
ts) 

tions (b-1) at p=13.5%. The length-scale parameters of the trials are chosen to match Pex P t( r )- 
The chord-distribution's (Table |j) indicate that model N (c=0) [shown in (b)] provides the 
best reconstruction. Each image has side length 39.7^m. 




5 10 . .15 20 



Figure 7: The auto-correlation function of the experimental image is well reproduced by the 
best reconstruction [model N (c=0)]. 



to the ice-water difference around 0°C, then a bulk modulus of 23.1GPa causes the N (c=0) 
model to agree perfectly with experiment at temperature points above the melting point of 
silver. 

The bounds are also shown in Fig. 11. For model N (c=0) the microstructure parameters 
are £i = 0.31 and 771 = 0.27. The results bound the experimental data and provide a reasonable 
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Figure 8: The good agreement between the model and experimental chord distributions 
[silver:/^ 1 ) (r); tungsten:^ 2 ) (r)] indicates that model N (c=0) provides the best reconstruc- 
tion. 



Figure 9: The best reconstruction [model N (c=0)]. The region shown is 198. 7x 198.7//m (cf. 
Fig. |). 



prediction of the Young's modulus below the melting point of silver. Note that even if the 
silver phase is given a non-zero bulk modulus past its melting point, the zero shear modulus 
causes the lower bounds for shear modulus and therefore Young's modulus to be identically 
zero. Unfortunately, there was no reported Poisson's ratio results for the composite, so we 
cannot compare to the model results for this quantity. 
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Figure 10: The silver phase (shown as solid) of the best reconstruction [model N (c=0)]. The 
side length is 39.7/mi. 
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Figure 11: The experimental Young's moduli compared with FEM data obtained from the best 
reconstruction [Fig. [1C[] . The Beran, Molyneux, Milton and Phan-Thien (BMMP) bounds for 
the reconstruction and Hashin and Shtrikman bounds are also shown. 



6 Analysis of analytical effective elastic moduli predictions 

We now compare the different analytical predictions of effective moduli with the finite element 
data. We have chosen to study the two most common models and only consider the moduli 



appropriate for the W-Ag composite studied above. The results are shown in Fig. 12(a) for 



overlapping spheres and in Fig. 12(b) for the single level-cut GRF or excursion set of Quiblier 
[model N(c=0)]. The self-consistent method provides a very good estimate of E e for model 
N(c=0), but not for overlapping spheres. This might be expected because the N(c=0) model 
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is symmetric with respect to phase-interchange (like the SCM) while the overlapping sphere 
model is not. As stated above the application of the generalized SCM is difficult because it 
is not obvious which phase should be chosen as the 'inclusion' phase. For the overlapping 
sphere model the tungsten phase is comprised of spheres (at 80% volume fraction), so is the 
more likely choice for the inclusion phase. Nevertheless, we report both estimates (80% W 
inclusions and 20% Ag matrix or 20% Ag inclusions and 80% W matrix) for both models. For 
either choice, the GSCM fails to provide an accurate estimate. Indeed, above the melting point 
of silver the GSCM vanishes for the case 20% Ag matrix case since the matrix phase is now 
completely soft. For the overlapping sphere model the Beran, Molyneux, Milton and Phan- 
Thien bounds are calculated using the microstructure parameters Ci = 0.52, r/i = 0.42 JH], 
Below the melting point of silver (where the contrast between the phases is moderate) the 
upper bounds provide a very good estimate of the effective moduli. 

A brief discussion of the effect of elastic contrast is necessary here. We have already noted 
that the analytical predictions of effective moduli do not explicitly depend on microstructure, 
but have a "built-in" microstructure. The elastic contrast, the ratio between the phase moduli, 
will determine how sensitive the effective moduli actually are to microstructure. For example, 
in the case of a two-phase composite having equal shear moduli but different bulk moduli, 
there is a simple exact formula for the effective bulk modulus which is totally insensitive 
to microstructure ||21||. In the case of small contrast, the effective moduli can be expressed 



exactly as a power series in the moduli differences [42]. Up to second order in this difference, 
at any volume fraction, the coefficients of the power series are not dependent on anything 
but the volume fractions and the individual phase properties. Therefore at small contrast, 
analytical predictions of effective moduli that explicitly depend only on volume fractions and 
phase moduli should all work well. 
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Figure 12: Comparison of theory (predictions and bounds) with finite element (FEM) calcula- 
tions for the Young's modulus of (a) the overlapping sphere model and (b) the single-cut GRF 
model [N(c=0), see Fig. [H|. The standard (SCM) and generalized (GSCM) self-consistent 
methods are shown, as are the Hashin and Shtrikman (HS) and Beran, Molyneux, Milton, 
and Phan-Thien (BMMP) bounds. 
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7 Conclusion 



Throughout this paper, we have treated the finite element computation method as being 
perfectly accurate, so that comparisons of elastic results to experiment were solely a test of 
how well the reconstructed microstructure compared to the real microstructure. This is not 
exactly true, since there are numerical errors in the finite element method [18, 17]. These 



are small, however, and are generally of about the same size or less than the differences seen 
between model computations and experimental data for the elastic moduli. There are also 
statistical sampling errors associated with the finite size (~40//m) of the models we employ 
to estimate the elastic properties. Since this is much greater than the correlation length of 
the samples (~5//m - see Fig. 0) we again assume these errors to be small. Therefore, the 
good agreement between model prediction and experimental data seen in this paper is good 
evidence that the model considered is indeed capturing the main aspects of the experimental 
microstructure. 

We have compared various theoretical results to finite element computations of the effec- 
tive Young's modulus E e for non-particulate media: a W-Ag composite and two model media 
(overlapping spheres and a single-cut Gaussian random field) . The generalized self-consistent 
method (derived for particulate composites) did not provide a good estimate of E e for the bi- 
continuous materials considered here. The standard self-consistent method provided a good 
estimate for the single-cut GRF and W-Ag composite. Since the method predicts zero moduli 
for porosity above 50% but the solid phase of the single-cut GRF remains connected up to 



porosities of around 90% J37| such agreement cannot be general. Upper bounds, calculated 
using three-point statistical correlation functions, provided a good prediction at low contrast 
{E\jEi ~ 6) for each composite. When one of the phases was completely soft the bounds 
lost predictive value. Therefore, for general composites, it is important to employ numerical 
computations of the effective moduli. For accurate numerical prediction of composite prop- 
erties it is important that a realistic model be used. Model-based statistical reconstruction, 
based on the Joshi-Quiblier-Adler approach, appears to be a viable route for microstructural 
simulation. However, it is important that the models underlying the procedure be capable of 
mimicking the composite microstructure. We have shown how several different models can be 
employed to find a useful reconstruction. 
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